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We study the growth of ahgned domains in nematic liquid crystals. Results are obtained solving 
the Beris-Edwards equations of motion using the lattice Boltzmann approach. Spatial anisotropy 
in the domain growth is shown to be a consequence of the flow induced by the changing order 
parameter field (backfiow). The generalization of the results to the growth of a cylindrical domain, 
which involves the dynamics of a defect ring, is discussed. 
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I. INTRODUCTION 



Liquid crystals [1] are an ideal material for the study 
of topological defects due to the complex textures they 
create which are easily visible to the naked eye. As topo- 
logical defects arise in many situations the observed phe- 
nomena in liquid crystals can be used to test theories in 
other areas of physics from cosmic strings [2] to vortices 
in superfluid helium [3] . Contrary to the assumption in- 
herent in most previous studies of defect dynamics in 
liquid crystals [4-6] , in a recent Letter [8] we found that 
backfiow, the coupling between the order parameter and 
the velocity fields, has a significant effect on the motion 
of defects. In particular, the defect speed can depend 
strongly on the topological strength in two dimensions 
and on the sense of rotation of the director about the 
core in three dimensions. 

These defects were free, in the sense that they were 
in an unbounded system. However, it is much easier to 
study liquid crystals experimentally in a confined system. 
A very straightforward example is the geometry used for 
display devices. In such a display the liquid crystal is 
sandwiched between two plates. As the optical and elec- 
trical responses of the liquid crystal are coupled, one can 
apply an electric field between the two plates and directly 
observe the behaviour. 

The operational state of many such devices, including 
pi-cells [9] , is topologically distinct from its state at zero 
voltage. Before the device can be used, the operational 
state must be nucleated and grow to fill the display. (A 
typical cross section through a domain wall separating 
an operational state and a zero voltage state is shown 
in Fig. 2(b).) These interfaces between topologically dis- 
tinct states can behave differently than nematic-isotropic 
interfaces studied by other groups [7] . Recent experimen- 
tal work on pi-cells [9] has shown an unusual anisotropy 
in the domain growth: one side of a domain grows faster 
than the opposite side. This would appear to be very sim- 
ilar to the anisotropy observed in our simulations of free 
defects [8]. However there are important differences in 
this system related to the wall tilt angle, which can dom- 



inate the defect-defect interaction energy, and the driving 
force of the electric field. As such, in order to unambigu- 
ously characterize the observed anisotropy we need to 
directly study the growth of these domain walls-which 
incorporate topological defects-in a confined geometry. 

In additional to the technological applications, simi- 
lar devices have been proposed as ideal experimental re- 
alizations of two-dimensional (2D) Ising models (in the 
plane defined by the walls of the device). Coarsening 
of reverse tilt domains in liquid crystal cells with het- 
erogeneous alignment layers has been shown to be con- 
sistent with predictions of the random-bond Ising model 
[10] thus providing experimental confirmation of theoret- 
ical predictions for domain growth under conditions of 
quenched random disorder. 

The dynamics of a liquid crystal medium is often mod- 
eled by using the Ericksen-Leslie-Parodi equations of mo- 
tion [1,20]. These equations describe the state of the 
liquid crystal in terms of a director field n which is re- 
lated to the orientation of the typically long, thin, rod- 
like molecules which make up the liquid crystalline mate- 
rial. The Ericksen-Leslie-Parodi equations are restricted 
to an uniaxial order parameter field of constant magni- 
tude. Thus they cannot model the dynamics of topo- 
logical defects where in the defect core the magnitude of 
order has a steep gradient and the order parameter field 
is biaxial. (However, they provide a good description of 
the bulk away from the defect core.) 

In order to describe the hydrodynamics of topological 
defects correctly, we use the more complex Beris-Edwards 
formulation of liquid crystal hydrodynamics [11]. The 
propensity to order, as well as the direction along which 
the system orders are conveniently described by a tensor 
order parameter Q [1]. The Beris-Edwards equations al- 
low for variations in the magnitude of the nematic order 
parameter as well as biaxiality present in defect cores. 
They model both defect dynamics and the coupling be- 
tween the velocity field and the motion of the order pa- 
rameter. We use a recent lattice Boltzmann algorithm 
[13] which has been shown to successfully model the full 
Beris-Edwards equations. 

Our aim in this paper is to study the growth of a do- 
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main of a ncmatic liquid crystal at the expense of a sec- 
ond domain with a different director orientation. Defects 
form at the walls between domains and their dynamics 
is vital in controlling the rate of growth. We find that 
a spatial anisotropy in domain growth can result from 
backflow and discuss how the wall speed varies with the 
material parameters of the liquid crystal, such as viscos- 
ity and elastic constants, the geometry and the surface 
properties of the confining cell and a external electric 
field. The Beris-Edwards equations of motion are pre- 
sented in Section II, and the results of the domain growth 
are described in Section III and Section IV. In Section V 
we discuss the relevance of our results to the experiments 
on pi-cells [9]. An outline of the numerical algorithm is 
given in the Appendix. 



II. THE HYDRODYNAMIC EQUATIONS OF 
MOTION 



We summarize the formulation of liquid crystal hy- 
drodynamics described by Boris and Edwards [11], ex- 
tended to include an electric field and surface potentials. 
Similar models have been examined by a number of re- 
searchers [12]. The continuum equations of motion are 
written in terms of a tensor order parameter Q which 
is related to the direction of individual molecules m by 
Qa/3 = {niarhj3 — ^Saj3) where the angular brackets denote 
a coarse-grained average. (Greek indices will be used to 
represent Cartesian components of vectors and tensors 
and the usual summation over repeated indices will be 
assumed.) Q is a traceless symmetric tensor. Its largest 
eigenvalue, |g, < g < 1, describes the magnitude of the 
order. 

We first write down a Landau-de Gennes free energy 
which describes the equilibrium properties of the liquid 
crystal [1,14] 

= dV {f bulk + f el + f field} + / dS{fsurf}- 

Jv Js (I 



(1) 



hulk describes the bulk free energy 

fbulk = -^(1 - ^)Qa/3 ^QafiQ/ijQ') 



(2) 



For 7 = 2.7 there is a first-order transition from the 
isotropic to the nematic phase. The minimum of fbulk 
describes a uniaxial nematic with an order parameter 
of the form Q^^jj = q{nan0 — |(5a/3) where q is zero in 
the isotropic phase and has a finite value in the nematic 
phase, and n is the director field. 

fei is the analogue of the Frank elastic free energy den- 
sity 



fel = ^{daQffjf + ^{dccQoci){d0Ql3'y) + 



L3 

-^QatsidaQiejidlsQje) 



(3) 



This can be easily mapped to give the Frank elastic con- 
stants Ki, K2 and K3 [11]. In particular, the "one elastic 
constant" approximation, Ki = K2 = K3 corresponds to 
Li > and L2 = L3 = 0. 

For a uniaxial nematic, the dielectric constant is 
anisotropic measured along or perpendicular to the di- 
rector. The relation between the electric displacement D 
and field E is of the form [1] 



D = e_LE-F(e|| -e_L)(n-E)n. 



(4) 



More generally, the dependence of the dielectric constant 
on the order parameter is described by 



where 



ea/3 — -^^aQa^ + Em^a/S 



2q 
2 1 



(5) 

(6) 
(7) 



giving results consistent with Eqn.(4) for the uniaxial ne- 
matic. The electric contribution to the thermodynamic 
potential f field is 



f field = jr> 



(8) 



The equation of motion for the nematic order param- 
eter is [11] 



{dt + u ■ V)Q - S(W, Q) = FH 



(9) 



where F is a collective rotational diffusion constant. The 
first term on the left-hand side of equation (9) is the ma- 
terial derivative describing the usual time dependence of 
a quantity advected by a fiuid with velocity u. This is 
generalized by a second term 

S(W, Q) = (^A + 0)(Q + 1/3) + (Q + I/3)(^A - Q) 
-2e(Q + I/3)Tr(QW) (10) 

where A = (W W^)/2 and f2 = (W - W^) /2 are the 
symmetric part and the anti-symmetric part respectively 
of the velocity gradient tensor Waj3 = dpUa- S(W, Q) 
appears in the equation of motion because the order pa- 
rameter distribution can be both rotated and stretched 
by flow gradients. This is a consequence of the rod-like 
geometry of the liquid crystal molecules. ^ is a constant 
which depends on the molecular details of a given liquid 
crystal. 

The term on the right-hand side of Eqn.(9) describes 
the relaxation of the order parameter towards the min- 
imum of the free energy. The molecular field H which 
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provides the driving motion is related to the derivative 
of the free energy by 

= Hbulk + Hel + Hfleld (11) 

where 

Hbulk = -Ao{l - |)Q + Ao7 (Q^ - (I/3)TrQ^) 

-Ao7QTrQ^ (12) 

(Hel)a/3 = Li{dj'^Qa/3) 

+ yap{dM^y (13) 

(Hfleld)a/3 = l^^-^"^!^ 1^-^7^)' (1^) 

and Saf) is the Kronecker delta. We work in a two- 
dimensional cross section, assuming that the order pa- 
rameter does not change in the perpendicular direction 
(although the director may point out of the simulation 
plane). In addition, the symmetry and zero trace of Q is 
exploited for simplification. 

At the surfaces of the device we assume a pinning po- 
tential 

fsurf = las{Qa0-Q%f- (15) 

We typically take Q° of the form 

Qlis = QiKnl - Sc^^m, (16) 

where q is set to the equilibrium bulk value. This cor- 
responds to specifying a preferred direction n" for the 
director at the surface. There can be other terms in the 
surface free energy [15] and a complete treatment of sur- 
face dynamics can be quite rich [16]. However, in this 
paper we will be operating in the strong pinning limit 
(as large) so that the only effect of the pinning potential 
is to furnish an almost fixed value of Qap at the surface 
(equal to Qq). In all cases studied here, the results are 
insensitive to the precise value of as, so long as it is large 
enough to be in the strong pinning limit. 
The fluid momentum obeys the continuity 

dtp + dapUa = (17) 

and the Navier-Stokes equation 

p{dt + Upd[j)Ua = dpTap + dfjCTaji 

+r]dp{(l - :idpPo)dyU.,Sap + daUp + dpUa) (18) 

where p is the fluid density and rj = pTf /3 is an isotropic 
viscosity (which is controlled by the simulation parameter 



Tf described in the Appendix). The form of this equa- 
tion is not dissimilar to that for a simple fluid. However 
the details of the stress tensor reflect the additional com- 
plications of liquid crystal hydrodynamics. There is a 
symmetric contribution 

1 5^ 

+ '^^{QaP + ^So,0)Q^eHje ' OpQ^, (19) 

and an antisymmetric contribution 

= Qa-yHjp — Ha-yQ^p. (20) 

These additional terms can be mapped onto the Ericksen- 
Leslie equations to give the Leslie coefficients [13]. The 
background pressure Pq is constant in the simulations to 
a very good approximation (±1%). 

The differential equations for order parameter field 
Eqn. (9) and the flow fleld Eqn. (18) are coupled. The 
velocity fleld and its derivatives appear in the equation 
of motion for the order parameter Eqn. (9). Unless the 
flow fleld is zero, u = 0, the dynamics given by Eqn. (9) 
are not relaxational and hydrodynamics play an impor- 
tant role. Conversely, the order parameter field affects 
the dynamics of the flow fleld through the stress tensors 
(19) and (20) which appear in the Navier-Stokes equation 
(18) and depend on Q and H. This back-action of the 
order parameter field on the flow field is usually referred 
to as backflow. To study these equations we use a lat- 
tice Boltzmann algorithm summarized in the Appendix. 
Other than when explicitly stated, the simulation param- 
eters are those listed in [17]. 



III. DOMAIN GROWTH 

We consider a liquid crystal conflned between two 
planes a distance apart. The director fleld may take 
topologically distinct states depending on the boundary 
conditions and applied voltage. In the simulations we set 
the boundary condition (Q" in Eq. (15)) so as to give 
a tilt angle —9p between the director and the y axis at 
X = and +6p at x = L^. At zero applied voltage these 
conditions result in a global minimum free energy state 
with a splayed director configuration, or horizontal (H) 
state as shown in Figure 1(a). At high voltages, typically 
on the order of 6V, the H state is no longer the global 
minimum, and a bend configuration (vertical state) is ob- 
tained like the one shown in Figure 1(b). At intermediate 
voltages, the vertical (V) state is more relaxed as shown 
in Figure 1(c). 

As the H and V states are topologically distinct, the 
transition from V to H requires nucleation of H domains 
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and the generation of defects. The problem we wiU inves- 
tigate is the growth (or shrinking) of the H state within 
the V state. In particular, we are interested in how 
hydrodynamics affects the speed of the domain walls. 
This is partly motivated by the observation in Ref. [9] 
that the domain growth in a liquid crystal device can be 
anisotropic and the speculation that this may be due to 
hydrodynamics. 

We have previously observed that the velocity of de- 
fects in unbound systems can be affected by hydrodynam- 
ics. In particular the defect speed can depend strongly 
on the topological strength in two dimensions and on the 
sense of rotation of the director about the core in three 
dimensions [8]. The crucial difference between the do- 
main growth problem and the motion of free defects is 
that in the latter case each defect moves due the direc- 
tor field of the other. In the domain growth problem 
the defects are not interacting but are dragged by the 
free-energy-driven movement of the domain walls. The 
free defects are accelerated as they approach each other 
while in the domain growth problem the defects move 
with a constant speed. Due to these differences and the 
additional geometrical parameters, a separate analysis is 
needed for the confined system which is also easier to re- 
alize experimentally and provides a better control of the 
parameters influencing the defect speed. 
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(a) (b) (c) 

FIG. 1. The possible alignment of directors for a tilt an- 
gle —6p on the top surface and +0p on the bottom surface 
{9p < 45 (leg: the surface tilt angle is measured with respect to 
the horizontal axis): (a) Director configuration when the field 
is switched off and the system had time to relax to its global 
minimum (H, or horizontal state); (b) the field is switched on 
at a fairly high voltage ~ 6V (V, or vertical state); (c) the 
field is at a voltage ~ 2V or lower. The system may remain 
in the metastable state (c) for some time even at zero volt- 
age. Periodic boundary conditions apply in the horizontal (y) 
direction. 

In order to study the role of hydrodynamics in the 
system we will examine the factors affecting the domain 
growth so that we can clearly identify what causes the 
wall speed anisotropy. The key parameters are the sur- 
face director tilt 6p, the sample thickness and material 



parameters such as coefficients in the bulk free energy (1): 
Aq, 7, and elastic constants Li, L2 and L3. In addition, 
the rotational diffusion constant F, which appears in the 
dynamical equation (9) for the order parameter gives an 
overall (inverse) time scale and is related to the Leslie- 
Ericksen viscosities [13]. 
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FIG. 2. (a) The initial director configuration is a horizon- 
tal domain (H state) in an otherwise vertically aligned sys- 
tem (V state), (b) As the system begins to relax, two defects 
are formed at the boundary of the horizontal and vertical 
domains. The left (right) defect has a topological strength 
s = —\ (s = -|-|). The curved arrows indicate the direction 
of the vortices induced by the reorientation of the director 
during the growth of the horizontal domain. The straight 
arrows point into the direction of defect motion. Note that 
there are periodic boundary conditions in the horizontal (y) 
direction. 

For simplicity, we will first study the undriven case 
of an H domain growing at the expense of a V state. 
We will later examine growth under the influence of an 
electric fleld. The initial conflguration, depicted in Fig. 
2(a), is a horizontal (i.e., along the y-direction) domain 
in an otherwise vertically aligned state. This models a 
time shortly after the electric field has been switched off 
when small but macroscopic domains have formed in the 
device. As the simulation proceeds, the director configu- 
ration relaxes rapidly to that shown in Fig. 2(b). During 
the relaxation defects are formed at the center of each do- 
main wall with strengths and — |, respectively. Once 
the two defects have formed the vertical domain begins 
to grow and the and —\ defects move in opposite 
directions. 

Our simulations correspond to a two-dimensional cross 
section of the two line defects, assuming that the order 
parameter does not change in the perpendicular direction 
(although the director may point out of the simulation 
plane). The two defects are topologically distinct only in 
two dimensions, but even in three dimensions they are 
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usually separated by an energy barrier. 

A particular advantage of the simulations is that the 
backflow can easily be switched off by setting = 
—PoSa/d and Tap to zero. (Compare this to (19) and 
(20).) Since there is no flow imposed, there is a zero 
velocity field throughout the whole simulation. The dy- 
namical equation in this case can be obtained from (9) 
by setting u to zero. It corresponds to the purely relax- 
ational Ginzburg-Landau model [18] 



dtQ = TH 



(21) 



where the molecular field H is given by (11). Compar- 
ing the dynamics obtained from the Ginzburg-Landau 
model and the full hydrodynamic equations, the effect of 
the backflow can be unambiguously identified. 

The Ginzburg-Landau equation (21) with a single elas- 
tic constant is invariant under a local coordinate trans- 
formation mirroring the director on the x axis. This cor- 
responds to the transformation 



Q 



xy 



Ixy, 



Q 



yx 



!yx- 



(22) 



The order parameter fields of the two moving defects with 
topological charges s = ±^ shown in Fig. 2(b) (even 
including the deformation due to the boundaries) trans- 
form into each other. Thus approaches based on a sim- 
ple Ginzburg-Landau equation predict that as the defects 
move they follow symmetric dynamical trajectories. 

We can construct a simple model for the domain mo- 
tion in the absence of hydrodynamic flow. In the bulk 
regions (away from the domain walls), fbuik is mini- 
mized with an uniaxial order parameter of the form 
Qap = Q.i'f^a'nfs — \5ai3)- We can then restrict our at- 
tention to the elastic free energy fej. With this form of 
the order parameter and if L2 — = 0, the elastic free 
energy density has the form fei = Liq^i^OY if the di- 
rector remains in the plane so that n = (cos 0, sin ^,0). 
The minima in the bulk regions (away from the domain 
walls) correspond to the director angle changing linearly 
along the x coordinate from —Op to +0p in the H domain, 
and from —Op to {+6p — 180 deg) in the V domain. The 
difference of the free energy densities of the two domains 
can then be written as 

A/ = /.-/H = ^(f-U (23) 

For Op < 45 deg, the horizontal domain grows because 

this decreases the free energy of the system. For Op > 
45 deg the horizontal domain should begin to shrink, and 
at Op = 45 deg the two domains have the same free energy 
and the defects should stop moving. 

If the H domain grows by a length of ALj, then the free 
energy of the system decreases by A/ x x A.Ly. Sim- 
ple relaxational arguments then suggest that the speed 
of domain growth can be described by the formula 
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X A/ X L^, 



(24) 



where rje is an effective viscosity. 

Surface tilt: We first investigate the effect of the sur- 
face tilt Op on the defect speed. Equations (23) and (24) 
suggest that as Op increases, the free energy difference 
decreases, and the defects should move more slowly. The 
defect speed v is plotted as a function of surface tilt Op in 
Fig. 3. Consider first the diamonds. These correspond 
to the case with backflow switched off. For this case both 
defects move at the same speed (but in opposite direc- 
tions). Notice that the defect velocity is proportional to 
{Op — 45 deg) . From (23) and (24) this leads to the con- 
clusion that the effective viscosity ije is independent of 
the tilt angle and its value is found to be 0.138 Pa s for 
the parameters of the simulation. 




FIG. 3. The velocity of the two defects as a function of sur- 
face tilt if backflow is ignored (diamonds) or included. Note 
that if backflow is not included then the two defects move 

with the same speed, which is well described by the dashed 
line based on Eqns. (23) and (24) Hydrodynamics accelerates 
the s = +^ defect (triangles) substantially, while it affects 
the s = — 5 defect (circles) much less. The speed anisotropy 
a is 36%. 

Back-flow: The triangles and circles in Fig. 3 show 
the velocity of the defects when backflow is included in 
the model. The s = +^ defect is considerably speeded 
up, whereas the = — 5 defect is only slightly ac- 
celerated. The defect speed reniains proportional to 
{Op — 45 deg) within a 2% error. The effective viscosities 
are ??+i/2 = 0.083 Pa s < r?_i/2 = 0.123 Pa s. These val- 
ues are comparable to the rotational viscosity 71 = 0.08 
Pa s [9]. The speed anisotropy defined as 

_ Au _ Vs=+l/2 - Vs=-l/2 
V {Vs=+l/2+Vs=-l/2)/^ 



(25) 



is independent of Op and its value is a = 36%. 

The order parameter field affects the flow field through 
the symmetric and antisymmetric stress tensors, (19) and 
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(20). The total nonviscous stress (i.e. the combination 
of all the stress terms not related to the velocity gradient 
tensor) is the sum of three terms 



cr + T = o-j 



CH + CTd- 



(26) 



Here <Ji,a/3 — —PoSaf3 is the stress due to the isotropic 
pressure. <7d,a/3 = — ^/^Q^^ ^g^Q is the deformation 
stress. For L2 = L3 = the deformation stress is 
Cd.a/J = —LiTr{daQdjsQ). The rest of the terms in (19) 
and (20) give what we will refer to as the molecular field 
stress, (Th, which is a function of H and Q. crd and the 
diagonal cTi do not change under the transformation (22) 
that transforms the defects of topological charge ±1/2 
into each other. Conversely, the off-diagonal elements of 
(Th have their sign inverted. Thus the stress fields and 
the resultant backflow are different for the two defects. 

The stress field era is related to the deformation free 
energy density, which is the same for both defects. It 
induces vortices similar to those around a solid cylinder 
moving in a viscous liquid. The flow points in the direc- 
tion of defect motion at the defect core. The contribution 
(Th describes the stress due to the the reorientation of 
the director. The reorientation is the strongest around 
the core while molecules near the surfaces reorient much 
less. The director reorientation induces vortices around 
the two defects as shown in Fig. 2(b). The direction of 
these vortices is determined by the gradient of the direc- 
tor angle taken moving around the defect in the positive 
direction. It is positive (negative) for the +^ (~^) de- 
fect. 
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(a) 



(b) 



FIG 

(b) s -- 



4. Velocity field corresponding to the (a) s = — ^ and 
= -1-^ defects shown in Fig. 2(b). There is a strong 



vortex pair around the 



+ i defect which, at the defect 



core, points in the direction of defect movement. The flow at 
the core of the s = — 5 defect is weaker and points opposite 
to the direction of defect propagation. 

The two contributions to the backflow reinforce each 
other for the .s = +i defect but partially cancel for the 
s = — i defect. The resulting flow fields can be seen in 
Fig. 4. The flow is stronger around the s = +h defect. 



Around the s = — ^ defect the flow is much weaker and 
the flow fleld points opposite to the defect propagation at 
the core. However, even in this case backflow accelerates 
the relaxational dynamics. 
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1.5 , 2 2.5 

(b) 

FIG. 5. (a) Speed of the -t-| (triangles) and — | (circles) 
defects as a function of the thickness of the sample. The di- 
amonds correspond to the case without hydrodynamics. The 
inset shows the relative speed anisotropy as a function of sam- 
ple thickness, (b) Effective viscosity r/e as a function of sample 
thickness for the case without hydrodynamics. The dashed 
lines in both figures correspond to the fit to the theoretical 
results discussed in the text. 



Scimple dimensions: As the sample becomes wider 
the speed of the defect propagation decreases [19] as can 
be seen in Fig. 5(a). The dependence of the defect veloc- 
ity on Lx follows from Eqns. (23) and (24) which give 
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V oc 



VeLx 



(27) 
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The effective viscosity can be calculated from ije = 
l/(27rsLir) J{V9q)'^dr, where 9q is the field due just to 
the defect itself [5,6] and the integral is over the volume of 
the system. Due to the confining geometry, one expects 
the integral to be dominated by the near-field contribu- 
tion (i.e. the field near the defect core) which is the same 
for both a static or moving defect [5]. The gradient of 9 
caused by a static defect drops off as 1/r and hence one 
expects the effective viscosity to go like \og{Lx/Lxo) [6]. 

This dependence is observed in Fig. 5(b). The fit is 
rje = 0.62 Pa s X log{Lx/LxQ) where Lr^o = 0.076/im is 
comparable to the defect core diameter. The dashed lino 
in Fig. 5(a) shows the fit to Eq. (27) including the 
dependence of the viscosity. 

When backflow is considered, the relative speed 
anisotropy increases with and saturates at about 60% 
as shown in the inset of Fig. 5(a). The increase is proba- 
bly due to the increasingly larger regions around the core 
involved in director reorientation. This leads to stronger 
vortices and hence stronger hydrodynamic effects. 
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FIG. 6. Speed of the +^ (triangles) and — ^ (circles) de- 
fects as a function of the external field V. Without hydrody- 
namics the defects move at the same speed (diamonds). The 
inset shows the relative speed anisotropy as a function of V. 

Electric Field: When an external electric field is ap- 
plied, it changes the free energy densities of the H and 
V domains. Thus it also influences the speed of the do- 
main growth. Fig. 6 shows the speed of the two defects 
as the function of applied voltage. For low voltages the 
free energy difference between the H and V domains can 
be estimated. If we assume that the orientation of the H 
and V states are unchanged from the zero voltage case 
(i.e. the director angle remains a linear function of x, as 
used in Eq. (23)), then substituting this into Equation 
(8) and integrating over space gives 



field 



1 



1 



ASttL^ \9p 7r/2 - 9j, 



V^LySm29p, (28) 



for a sample of length Ly. The total free energy differ- 
ence between the domains is the sum of the elastic and 
the field contributions, Eqn. (23) and Eqn. (28), re- 
spectively. Substituting this into Eqn. (24) results in a 
parabolic dependence of the velocity on V for low volt- 
ages as shown in Fig. 6. At Vumit ~ 0-9V domain growth 
is reversed since the effect of the surfaces is balanced by 
the influence of the electric fleld. 
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FIG. 7. A V domain growing in an asymmetric Ha envi- 
ronment. At high fields the horizontal domain is deformed 
moving the defects towards one of the boundaries. 

At high voltages (F ~ O.GF), the H domain is replaced 
by Ha, an asymmetric horizontal domain shown in Fig. 
7. It has a lower free energy than the H domain due to 
the more favorable alignment with the electric field. As 
a result of the deformation of the Ha domain, the defects 
move towards the surfaces. For V > W the molecules in 
the bulk are almost completely aligned with the vertical 
field. The horizontal region is then confined to a thin 
layer near the surfaces. For this type of configuration, 
the free energy difference between the domains is pro- 
portional to the voltage [9]. This gives the linear slope 
of the curve in Fig. 6 seen for the higher voltages. 

For low voltages the speed anisotropy is 36% and inde- 
pendent of the voltage. As shown in the inset of Fig. 6, at 
high voltages {V > IV) the speed anisotropy decreases. 
The effect on the anisotropy can be explained by the rel- 
ative weight of the relaxational dynamics and the back- 
fiow. In the ^ = 1^ to 1.6V range the relaxational 
dynamics are substantially speeded up by the increasing 
voltage due to the electric field contribution (8) in the 
free energy. The backflow is induced by the the stress 
fields given in Eq. (19) and Eq. (20). These stress fields 
do not depend directly on the electric field, only on the 
order parameter field, which changes only slightly in this 
voltage range. Therefore the stress fields do not increase 
with the increasing voltage, and the hydrodynamics is 
dominated by the relaxational dynamics at high fields. 
(In comparison, for the experiment presented in [9] the 
domain wall speed was ~ 0.1/um/ms and the anisotropy 
also decreased with increasing voltage.) 



IV. OTHER CONTROL PARAMETERS 

The equations governing the dynamics of liquid crys- 
tals covered in Section II contain a large set of parame- 
ters. In this section we explore some of this parameter 
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space. In particular, wc examine the case of multiple 
elastic constants, important for comparing to any real 
liquid crystal. In addition we look at the influence of the 
different viscosities and free energy parameters on the 
balance of the relaxational dynamics and the backflow. 
We also examine the case of asymmetric and inhomoge- 
neous surface tilts since they give new insights about the 
underlying symmetries of the system and are important 
for practical devices with non-trivial surface anchoring. 

Elastic constants: Real liquid crystals have multiple 
elastic constants. We now consider the effect of non-zero 
elastic constants L2 and L3. If L2 ^ in the expression 
of the molecular field Eqn. (13), the dynamical equa- 
tion in the absence of backflow Eq. (21) loses its invari- 
ancc under the transformation (22). However, the speed 
anisotropy is very small. The reason for this is that the 
relaxational dynamics are still invariant under the mir- 
roring transformation for a uniaxial order parameter with 
a constant magnitude [22] . In our setup these conditions 
hold except for close to the defect cores. 



0.45 




0.08 
r| (Pa s) 



0.14 



FIG. 8. Speed of the -1-^ (triangles) and — ^ (circles) de- 
fects as a function of the viscosity 77 = The dashed line 
corresponds to the velocity without backflow. The inset shows 
the relative speed anisotropy as a function of r). 



A larger anisotropy in speed is obtained if L3 7^ 0. If 
L3 > (La < 0) then the s = +^ {s = -i) defect is 
faster. For Li =8.73pN and L3 =15.87pN, we measure 
a speed anisotropy of a = 3%, in a system without hy- 
drodynamics. This anisotropy due to the unequal elastic 
constants increases with increasing sample thickness. For 
Lx = 1.25/im (our benchmark system has = 0.7/Ltm 
[17]) the anisotropy due to non-zero Z/3 is a = 6%. This 
may be due to the fact that, for free defects, the elas- 
tic anisotropy causes significant deviations from the case 
of isotropic elasticity in the order parameter field away 
from the axis determined by the two defect cores [1]. In 
the thinner sample, the surfaces "cut off' this part of the 
field, decreasing the anisotropy. 




FIG. 9. Speed of the -|-| (triangles) and — i (circles) de- 
fects as a function of F. The diamonds correspond to the case 
without hydrodynamics when v oc F. The inset shows the 
relative speed anisotropy as a function of F. 

If the surface tilt is close to vertical, and the horizon- 
tal domain is shrinking, then the L3 dependence of the 

speed anisotropy is the opposite. If < (L3 > 0) then 
the s = +i(s = — i) defect is faster. Since in the two 
cases (growing vs. shrinking domain) the order parame- 
ter fields near the axis of the two cores are the same, this 
should also be attributed to the differences in the order 
parameter fields far from the defect cores which results 
from the different tilts. 

Viscosities and diffusion: Consider now the effect 
of the parameters governing the time scales in the equa- 
tion of motion for the domain growth, r/ is proportional 
to the viscosity in the Navier-Stokes equation (18). In- 
creasing Tf increases the viscosity, slows the defects, and 
decreases their velocity anisotropy as shown in Fig. 8. 
The velocity tends to that measured without backflow, 
as represented by a dashed line in the figure. This is 
as expected since backfiow will be suppressed by a large 
viscosity. 

Increasing T which appears in Eq. (9) increases the 
velocity of both defects, but decreases the relative speed 
anisotropy as shown Fig. 9. The defects move faster, due 
to the fact that F governs the speed of the relaxational 
dynamics. Since the stress fields in Eq. (19) and Eq. (20) 
do not depend on F, they do not increase. As a result, as 
F increases, the relaxational dynamics speed up, but the 
backfiow does not change as significantly, and as a result 
the anisotropy decreases. 

Free energy: 7 is the parameter in the free energy 
(2) which controls the magnitude of order in the bulk 
of the domains. (The isotropic-nematic phase transition 
is at 7 = 2.7.) When 7 decreases, the defect core gets 
bigger and this results in a smaller effective viscosity [5] . 
Thus the defects move faster under the relaxational dy- 
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namics. The decrease in the magnitude of order results 
in a smaller backflow due to the reorientation and a de- 
creases. 
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(a) 




(b) (C) 

FIG. 10. (a) Director and velocity field about the (b) 
+1/2 and (c) —1/2 defects for asymmetric surface anchor- 
ing. The director tilt is Op{x = 0) = — lOdeg at the top and 
6p{x = Lx) = -1-60 deg at the bottom surface. The qualitative 
differences between the flow-fields of the two defects are the 
same as for the symmetric case in Fig. 4. 

Changing Aq in the free energy (2) does not affect the 
homogeneously aligned bulk H and V domains, only the 
defect structure. The larger Aq, the larger the energy 
cost of any deviation from the magnitude of order corre- 
sponding to the minimum of fbuik in (2). Decreasing Aq 
increases the size of the defects, and as above, the defects 
move faster under the relaxational dynamics. Since the 
defect core size increases, the magnitude of order around 
the core decreases resulting a smaller backflow due to re- 
orientation and hence a decreases. The effect of increas- 
ing Li is similar. It increases the defect core size, speeds 
up the defects, and leads to a smaller velocity anisotropy. 

Non-symmetric surface tilt: The director surface 
tilt at the top and bottom surfaces does not have to be 
symmetric. In this more general case, tilt angles at the 
surfaces can be written as a sum of a symmetric and an 
antisymmetric contribution: 



9p{x = 0) 
9p{x = Lx) 



+ 0: 



p,s 



(29) 



The dynamical equations for L2 = = and with- 
out flow are invariant under the local rotation of all the 



molecules by the same angle. The dynamics for a non- 
symmetric surface tilt can therefore be obtained from 
the symmetric case by rotating all the molecules by Op^a- 
Thus even for non-symmetric tilts the two defects move 
with the same speed and remain in the center between 
the two plates. The full dynamical equations however, 
are not invariant under the rotation of molecules if flow 
is included (or if L2 7^ 0, L3 7^ 0). Thus the defects 
will move with a different speed. Moreover they are no 
longer constrained by symmetry to lie midway between 
the two plates. Typical director and velocity flelds for an 
asymmetric surface tilt arc shown in Fig. 10. 

It is also possible to construct a patterned alignment 
of liquid crystals on surfaces [23]. If the surface tilt is 
not homogeneous then, when the defect arrives at a re- 
gion with a different tilt, it assumes the new velocity 
corresponding to this tilt. On the boundary of two bend 
domains with opposite surface tilt the defect can even 
"change" topological strength by merging with the two 
defects located at the surface as shown in Fig. 11. 
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(b) 

FIG. 11. The influence of surface tilt inhomogeneity; the 
surface tilt changes from 6p = +15dcg to Op = —15 deg to- 
wards the right hand side of the figure, (a) Upon reaching the 
border of the two bend domains with opposite surface tilt, the 
defect +i merges with the the two — i defects located at the 
surfaces; (b) a — ^ defect is formed. 



V. GROWTH OF A CYLINDRICAL DOMAIN: 
THE HYDRODYNAMICS OF A DEFECT RING 

In this section we consider the three dimensional ana- 
logue of our system, where a liquid crystal is held be- 
tween parallel plates iim apart. A domain nucleated 
at a point will grow to a cylindrical shape with its axis 
perpendicular to the confining plates. Based on our ex- 
perience with the two dimensional system wc can discuss 
how backflow can affect the growth of a cylindrical H 
domain in a V environment. 
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At the domain boundary there is a defect ring, as 
shown in Fig. 12(a). The defect configuration of a verti- 
cal cross section through the ring (perpendicular to the 
plates) changes gradually from a —5 to a +5 defect. Al- 
though this problem is three-dimensional, a vertical cross 
section indicated by dashed lines in Fig. 12(a) gives a ge- 
ometry similar to that considered in Section IV. 




Rubbing direction 

(a) 




FIG. 12. (a) Confined between two horizontal surfaces, a 
cylindrical H domain is growing in a V environment. There is 
a defect ring (dotted line) at the domain boundary, (b) The 
cross sections indicated by dashed lines in (a) are shown. This 
corresponds to the simulation plane considered in this paper. 
The director and the vortices due to the reorientation axe in 
the plane of the cross section, (c) The cross section indicated 
by dashed-dotted lines in (a). The directors are perpendicular 
to the plane of the cross section. The plane of vortices due 
to the reorientation is also perpendicular to the cross section. 
However, for both (b) and (c) the vortices due to the defor- 
mation stress are in the plane of the cross section and the 
flow at the core points in the direction of defect propagation 
indicated by straight arrows. 

For simplicity, assume that the H domain is a perfect 
cylinder. In this case the director field of any vertical 
cross section passing through the middle of this cylin- 
drical domain can be obtained from our two-dimensional 
simulation plane by rotating the director field locally by 
a given angle around the vertical axis, x. 

Let us now examine the effect of the backflow. cr^ and 
the diagonal ai do not change during the local rotation 
of the tensor order parameter field Q around the vertical 
axis, X, by the same angle. The generated velocity vor- 
tices always lie in the plane of the cross sections across 



the defects ring. These generate flow which is always in 
such a direction as to expand the defect ring, indicated 
by straight arrows in Figs. 12(a-c). 

The tensor cth does however change under a local ro- 
tation around x. For the cross section which includes 
a -|-^ defect the resultant backflow points in the direc- 
tion of defect motion. For the —\ defect, the flow points 
opposite to the defect motion, as shown in Fig. 12(b). 
In both cases the vortices are in the plane of the cross 
section. Fig. 12(c) shows the cross section indicated by 
the dashed-dotted lines in Fig. 12(a). Now the directors 
and the vortices due to the reorientation are in a plane 
perpendicular to the cross section. Thus the total flow 
field will vary around the domain wall and will lead to 
anisotropy in the domain growth. 

An experimental setup similar to this was consid- 
ered by Acosta et al. in their investigation of domain 
growth and switching in pi-cell liquid crystal devices [9]. 
The growth of a horizontal domain in a bend (V) or 
twisted bend environment was studied: such a transition 
is needed to produce the operational state of the device. 
(The twisted bend configuration has a lower energy than 
the bend state for small surface tilts and low voltage, if 
the Frank elastic constant K2 is sufficiently small. The 
twisted bend state is replaced by a bend state for larger 
(^ 30deg) surface tilts.) A cylindrical bend or twisted 
bend domain was formed in a H environment and the do- 
main wall velocity was measured at four points around 
the ring, where its cross section corresponds to a +i and 
a — ^ defect, and at two points halfway between these. 
It was found that the wall at the +\ defect moved sub- 
stantially faster than at the other three. It seems very 
plausible that the essential physics is captured by our 
model. 

Further measurements of defect dynamics in confined 
geometries have been done very recently [24] and these 
techniques should allow further testing of the concepts 
we present here. 



VI. SUMMARY 

In this paper we explored domain growth in nematic 
liquid crystals. Defects form at moving domain walls. 
We find that a wall incorporating a s = +^ defect is 
substantially speeded up by backflow effects, whereas a 
wall containing a s = — 5 defect is only slightly affected. 
This was explained in terms of the symmetry properties 
of the different stresses acting on the defects. These re- 
inforce each other for the s = +^ defect while partially 
cancelling for the s = —\ defect. The influence of dif- 
ferent material and geometrical parameters on the veloc- 
ity anisotropy was interpreted by comparing the relative 
weight of the relaxational dynamics and the backflow. By 
generalizing two-dimensional simulation results, a quali- 
tative picture was proposed for the role of the backflow 
in three dimensions. 
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Results were obtained using a lattice Boltzmann al- 
gorithm to solve the Beris-Edwards equations of liquid 
crystal hydrodynamics. Working within the framework 
of a variable tensor order parameter it was possible to 
correctly incorporate variations in the magnitude of or- 
der and hence the dynamics of domain walls and their 
associated topological defects. 
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APPENDIX A: A LATTICE BOLTZMANN 
ALGORITHM FOR LIQUID CRYSTAL 
HYDRODYNAMICS 

We now summarize a lattice Boltzmann algorithm 
which solves the hydrodynamic equations of motion of 
a liquid crystal (9), (17), and (18). Lattice Boltzmann 
algorithms are defined in terms of a set of continuous 
variables, usefully termed partial distribution functions, 
which move on a lattice in discrete space and time [21]. 

The simplest lattice Boltzmann algorithm, which de- 
scribes the Navier-Stokes equations of a simple fluid, is 
defined in terms of a single set of partial distribution 
functions which sum on each site to give the density. For 
liquid crystal hydrodynamics this must be supplemented 
by a second set, which are tensor variables, and which 
are related to the tensor order parameter Q [13]. 

Wc define two distribution functions, the scalars fi{x) 
and the symmetric traceless tensors {x) on each lattice 
site X. Each /j, Gi is associated with a lattice vector ej. 
We choose a nine-velocity model on a square lattice with 
velocity vectors = (±1, 0), (0, ±1), (±1, ±1), (0, 0). 
Physical variables are defined as moments of the distri- 
bution function 

i i i (Al) 

The distribution functions evolve in a time step At 
according to 

fi{x + eiM,t + M)- fi{x,t) = 

^ [Cfi{x, t, {/O) + Cfi{x + eAt, t + At, {/*})] , 

(A2) 

Gi{x + eiAt,t + At) - Gi{x,t) = 
At 

— [Cciix, t, {Gi}) + Coiix + eiAt, t + At, {G*})] . 

(A3) 



This represents free streaming with velocity and a col- 
lision step which allows the distribution to relax towards 
equilibrium. /* and G* are first order approximations to 
fi {x + Si At, t + At) and Gj {x + CiAt, t + At) respectively. 
They are obtained by using only the collision operator 
Cfi{x, t, {fi}) on the right of Equation (A2) and a simi- 
lar substitution in (A3). Discretizing in this way, which 
is similar to a predictor-corrector scheme, has the advan- 
tages that lattice viscosity terms are eliminated to second 
order and that the stability of the scheme is improved. 

The collision operators are taken to have the form of a 
single relaxation time Boltzmann equation [21], together 
with a forcing term 

cMx,t,{fi}) = --iMx,t) - f!\x,t,{f,})) 

+p,{x,t,{f,}), (A4) 
CGi{x,t,{Gi}) = --{Gi{x,t) - Gf (f,t,{Gi})) 

+Mi{x,t,{Gi}). (A5) 

The form of the equations of motion and thermodynamic 
equilibrium follow from the choice of the moments of the 
equilibrium distributions /f* and G^' and the driving 
terms pi and Mj. Full details of the algorithm can be 
found in [13]. 
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